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Abstract:  Despite  the  critical  importance  of  resolving  flow  at  or  near  domain  boundaries,  often 
boundary  condition  formulations  are  implemented  in  an  ad-hoc  fashion.  The  purpose  of  this  work 
is  to  provide  a  general  framework  for  the  implementation  and,  importantly,  the  verification  of 
boundary  conditions  for  node-  and  cell-centered  finite  volume  schemes.  These  conditions  may 
include  any  combination  of  Dirichlet  conditions,  Neumann  conditions,  extrapolation,  and  in  some 
instances,  the  conservation  equations  themselves.  Specific  conditions  for  inviscid  walls,  inflow,  and 
outflow  are  systematically  tested  using  manufactured  and  exact  solutions.  The  procedures  in  this 
work  provide  a  method  for  the  verification  of  complex  boundary  conditions  as  well  as  shed  valuable 
physical  insight  for  different  problem  situations. 

Keywords:  Boundary  Conditions,  Finite  Volume,  Lagrange  Multipliers,  Manufactured  Solutions, 
Verification. 


1  Introduction 

Computational  fluid  dynamics  (CFD)  is  used  to  simulate  flows  for  a  variety  of  applications.  Often,  these 
applications  require  high  levels  of  accuracy  at  or  near  domain  boundaries.  Examples  of  such  applications 
include  calculations  of  aerodynamic  lift  and  drag,  computation  of  boundary  layer  profiles,  estimation  of 
heating  at  surface  locations,  and  mass  flux  computation.  For  these  and  other  cases,  the  region  near  the 
boundary  is  the  primary  focus  of  the  CFD  simulation.  Although  the  computational  domain  may  extend 
greatly  from  these  boundaries,  the  areas  near  them  often  require  the  greatest  resolution  and  numerical 
accuracy. 

Despite  the  need  for  high  accuracy  near  boundaries  (or  perhaps  because  of  it),  numerous  boundary 
procedures  have  been  proposed  in  the  literature.  For  example,  a  widely  used  inviscid  wall  treatment  by 
Rizzi  [1]  involves  the  use  of  the  momentum  equation  to  obtain  the  pressure.  Jameson  directly  modified 
the  flux  at  an  inviscid  wall  to  involve  no  convective  contribution  [2].  Alternately,  Dadone  and  Grossman 
advocate  a  curvature  corrected  symmetry  condition  for  an  inviscid  wall  [3].  Balakrishnan  and  Fernandez 
advocate  a  variety  of  other  methods  for  an  inviscid  wall  involving  additional  quantities  such  as  entropy  and 
enthalpy  [4].  Numerous  other  strategies  exist  for  inviscid  walls  as  well  as  other  boundary  conditions,  such 
as  inflow,  outflow,  and  no-slip  walls.  The  difficulty  is  that  many  of  these  methods  have  not  been  rigorously 
verified  and  may  or  may  not  be  consistent  with  the  interior  discretization  schemes.  Recently,  Choudhary  et. 
al.  [5]  successfully  verified  boundary  conditions  using  the  method  of  manufactured  solutions  (MMS).  Their 
approach  requires  carefully  constructed  manufactured  solutions  that  already  satisfy  the  boundary  conditions, 
which  means  that  a  new  manufactured  solution  needs  to  be  constructed  for  each  boundary  condition  and 
geometry. 

This  work  addresses  many  of  these  difficulties  by  focusing  on  three  main  goals.  First,  we  provide  a 
general  framework  for  implementation  of  boundary  conditions  for  both  node-  and  cell-centered  schemes.  The 
framework  applies  to  arbitrary  boundary  conditions  for  any  physical  system.  Second,  we  provide  a  rigorous 
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boundary  nodes  (node-centered) 

(a)  Boundary  degrees  of  freedom  for  node-centered  methods. 


ghost  nodes  (cell-centered) 

(b)  Boundary  ghost  nodes  for  cell-centered  methods. 


Figure  1:  Location  of  boundary  condition  enforcement  for  node-  and  cell-centered  schemes. 


and  general  approach  for  the  verification  of  these  boundary  conditions  based  on  MMS.  While  verification  is 
an  essential  step  in  CFD  development  and  practice,  many  verification  strategies  have  neglected  boundary 
effects  [6,  7,  8,  9,  10].  In  this  work  we  extend  the  traditional  use  of  MMS  to  include  boundary  conditions 
and  uncover  certain  key  aspects  of  using  MMS  in  this  way.  While  the  verification  procedure  determines  the 
numerical  accuracy  of  the  discretization,  it  does  not  indicate  that  the  proper  equations  have  been  selected  for 
a  given  problem.  Therefore,  the  third  objective  of  this  work  is  to  demonstrate  the  requirement  of  choosing 
physically  correct  boundary  conditions  for  a  given  problem.  Unlike  the  interior,  at  boundaries,  the  choice  of 
governing  equations  is  not  well  established.  The  studies  in  this  work  indicate  that  the  choice  of  boundary 
conditions  is  not  unique,  but  careful  analysis  is  necessary  to  verify  the  properties  of  these  formulations. 

In  order  to  accomplish  these  tasks,  we  explore  the  formulation  of  a  variety  of  common  boundary  types, 
including  inviscid  walls,  inflow,  and  outflow,  for  both  node-  and  cell-centered  finite  volume  schemes.  For  each 
of  these  boundary  types,  we  explore  a  variety  of  conditions  spanning  Dirichlet,  Neumann,  and  extrapolation. 
Importantly,  we  highlight  the  fact  that  node-centered  schemes  allow  for  direct  use  of  the  governing  equations 
of  motion  (mass,  momentum,  and  energy)  themselves  at  the  boundaries.  Cell-centered  schemes  do  not  easily 
allow  use  of  the  equations  of  motion,  but  usually  rely  on  extrapolation  or  other  conditions.  Great  care 
must  be  taken  with  cell-centered  schemes  to  avoid  the  selection  of  additional  conditions,  which  may  conflict 
with  the  interior  equations  of  motion.  Other  boundary  condition  types  and  implementations  beyond  those 
discussed  in  this  work  are  certainly  possible.  Here,  the  principal  objective  is  to  provide  a  framework  for  the 
implementation  and  verification  of  any  boundary  condition. 

The  paper  is  organized  as  follows:  First,  boundary  condition  formulations  are  explored  for  node-  and 
cell-centered  finite  volume  schemes.  We  present  a  common  discretization  framework  with  which  we  test  a 
variety  of  governing  boundary  equations.  We  show  how  to  verify  these  boundary  conditions  using  MMS. 
We  then  present  a  number  of  grid  refinement  studies.  First  a  quasi-ID  nozzle  is  tested  with  both  MMS  and 
exact  solutions.  Next,  we  extend  the  use  of  MMS  to  two  dimensions  and  wall  boundary  conditions.  Finally 
we  present  error  convergence  results  using  Ringleb  flow  and  a  NACA  0012  subsonic  airfoil.  We  then  offer 
some  conclusions  derived  from  these  studies. 


2  Boundary  Condition  Implementation 


In  this  work  we  test  a  variety  of  boundary  condition  implementations  and  assess  the  resulting  impact  on 
accuracy  through  rigorous  verification  studies.  In  the  interior  of  the  domain,  we  solve  the  steady  Euler 
equations, 


§?  +  v 

OT 


F  =  0. 


(1) 


where  Q  =  (p,  pu,  pv,  pe)T  is  the  vector  of  conserved  variables,  and  F  is  the  inviscid  flux  vector.  Here,  p  is 
density,  u  and  v  are  the  Cartesian  velocity  components,  and  e  is  the  total  energy  per  unit  mass.  We  are 
interested  in  solving  the  steady  equations  to  which  we  add  the  pseudo-time  (r)  derivative  for  convenience  in 
marching  to  steady  state.  We  test  both  node-  and  cell-centered  finite  volume  spatial  discretizations,  which 
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result  in  a  semi-discrete  set  of  non-linear  equations  of  the  form 


+  R(Q)  =  0,  (2) 

where  R(Q)  is  the  steady  residual  at  either  an  interior  node  or  an  interior  cell  location  depending  on  the 
discretization  procedure. 

In  addition  to  the  interior  discretization  scheme,  we  must  also  incorporate  boundary  conditions  to  close  the 
system  of  equations  on  a  finite  domain.  Here  we  focus  on  inviscid  wall,  inflow,  and  outflow  conditions  in  order 
to  explore  fundamental  issues  of  physical  and  numerical  accuracy.  The  methodology  developed  here  is  quite 
general  and  directly  applies  to  other  boundary  condition  types  as  well.  For  node-centered  discretizations, 
boundary  conditions  are  enforced  directly  at  the  boundary  nodes  coincident  with  the  physical  boundary, 
shown  in  Figure  1(a).  For  cell-centered  discretizations,  we  introduce  additional  unknowns  in  the  form  of 
ghost  nodes  located  at  the  flux  quadrature  points  of  the  boundary  faces,  shown  in  Figure  1(b).  The  ghost 
nodes  are  then  used  in  an  upwind  flux  formula  to  determine  the  numerical  flux  through  the  boundary  face.  In 
this  manner,  the  cell-centered  boundary  formulation  remains  water-tight.  The  node-centered  configuration, 
however,  is  not  strictly  water-tight  since  the  fluxes  surrounding  the  boundary  nodes  do  not  always  cancel 
with  nearby  interior  nodes. 

For  both  node-  and  cell-centered  formulations,  we  define  a  “boundary  residual,”  Rb(Q),  which  we  drive 
to  zero  at  steady  state  along  with  the  interior  residuals,  R(Q): 

Rb(Q )  =  o.  (3) 

At  boundaries  in  a  node-centered  discretization,  replaces  R  at  the  boundary  nodes.  In  cell-centered 
discretizations,  Rb  provides  the  governing  equations  for  the  ghost  nodes.  In  all,  we  test  fifteen  different 
boundary  implementations,  which  are  listed  with  a  common  notational  convention  in  Table  1  for  clarity.  The 
methods  all  involve  Dirichlet-specified  quantities,  with  the  state  specification  completed  by  other  methods, 
such  as  Neumann  conditions,  extrapolation,  or  combinations  of  the  equations  of  motion.  We  will  refer  to 
this  table  as  we  develop  various  forms  for  Rb  in  the  following  sections. 


Table  1:  Notation  for  boundary  condition  methods  tested. 


Short  Name 

Boundary  Type 

Dirichlet 

Other  conditions 

n-INVl 

Node-centered 

Inviscid  wall 

Un 

Lagrange  multipliers 

n-INV2 

Node-centered 

Inviscid  wall 

Un 

Ut  Neumann,  mass  and  energy  eqs. 

n-INV3 

Node-centered 

Inviscid  wall 

Un 

p,  Ut,  pe  Neumann 

n-INV4 

Node-centered 

Inviscid  wall 

- 

zero  convective  flux 

n-INFl 

Node-centered 

Inflow 

s,h,ut 

Lagrange  multipliers 

n-INF2 

Node-centered 

Inflow 

s,h,ut 

outgoing  characteristic  eq. 

n-OUTl 

Node-centered 

Outflow 

P 

Lagrange  multipliers 

n-OUT2 

Node-centered 

Outflow 

P 

outgoing  characteristic  eqs. 

c-INVl 

Cell-centered 

Inviscid  wall 

p ,  ut,  pe  extrapolated 

C-INV2 

Cell-centered 

Inviscid  wall 

Pressure  extrapolated 

C-INV3 

Cell-centered 

Inviscid  wall 

p,ut,  pe  Neumann 

c-INFl 

Cell-centered 

Inflow 

s,  h,  ut 

un  extrapolated 

C-INF2 

Cell-centered 

Inflow 

R~,ut,s 

R+  extrapolated 

c-OUTl 

Cell-centered 

Outflow 

P 

p,  u,  v  extrapolated 

C-OUT2 

Cell-centered 

Outflow 

R~ 

R+  ,Ut,  s  extrapolated 

2.1  Node-Centered  Boundaries 

All  boundary  conditions  involve  the  specification  of  a  certain  number  of  Dirichlet  (or  Neumann)  conditions 
augmented  by  additional  information  derived  from  the  interior  field.  With  node-centered  schemes  it  is 
straightforward  to  select  a  combination  of  the  governing  equations  of  motion  (mass,  momentum,  and  energy) 
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to  enforce  at  boundary  nodes.  This  is  because  the  boundary  nodes  lie  within  control  volumes  for  which  we 
can  easily  implement  flux  balances.  This  will  be  shown  in  the  methods  described  below.  In  contrast,  cell- 
centered  boundary  conditions  often  require  the  use  of  ghost  nodes  for  which  there  is  no  natural  control 
volume  or  flux  balance.  Thus,  it  becomes  difficult  to  apply  the  equations  of  motion  directly  at  the  ghost 
nodes.  Instead,  we  choose  some  other  method  to  define  the  ghost  node  state,  such  as  extrapolation. 

Enforcement  of  boundary  conditions  within  a  node-centered  scheme  involves  the  specification  of  the 
boundary  residual,  J4,  directly  at  the  boundary  nodes  shown  in  Figure  1(a).  In  this  section,  we  outline  two 
procedures  to  obtain  the  boundary  residual.  The  first  involves  the  use  of  Lagrange  multipliers  as  discussed 
in  detail  by  Allmaras  [11].  The  second  involves  the  multiplication  of  the  governing  equations  of  mass, 
momentum,  and  energy  by  a  selection  matrix  to  complete  the  boundary  conditions.  In  addition,  we  discuss 
a  third  method  involving  a  commonly  used  weak  boundary  condition  for  an  inviscid  wall. 


2.1.1  Lagrange  Multipliers 

A  method  of  boundary  condition  enforcement  that  has  enjoyed  widespread  use  in  the  finite  element  commu¬ 
nity  for  several  decades  is  the  Lagrange  multiplier  method  introduced  by  Babuska  [12].  The  extension  of  this 
method  to  the  Navier-Stokes  equations  was  discussed  in  detail  by  Allmaras  [11],  but  seems  to  have  received 
little  attention  in  the  finite  volume  community.  The  reader  is  reffered  to  the  above  sources  for  details.  It 
suffices  to  state  here  that  the  method  involves  a  modification  of  the  variational  statement  to  include  extra 
conditions  along  Dirichlet  boundaries,  leading  to  an  extended  system  of  the  form 


(4) 


Here,  B(Q)  —  b  represents  a  vector  of  to  additional  Dirichlet  conditions  that  must  be  satisfied  at  boundary 
nodes.  For  example,  at  a  fixed  inviscid  wall,  m  =  1,  B(Q)  =  u  ■  n,  and  6  =  0.  The  additional  m  conditions 
are  incorporated  via  a  vector  of  m  Lagrange  multipliers,  A.  The  choice  of  variables  for  Q  in  Equation  4 
produces  different  conditions.  While  Allmaras  chooses  entropy  variables,  Q  =  —  e,u,v,  —  1)T, 

other  variables,  such  as  the  primitive  variables,  Q  =  (P,u,v,T)t ,  are  possible,  leading  to  different  boundary 
equations.  Thus,  the  Lagrange  multiplier  method  does  not  appear  to  uniquely  define  the  boundary  equations. 

While  the  extended  system  in  Equation  4  can  be  solved  directly  to  include  the  Lagrange  multipliers, 
these  may  be  eliminated  from  the  problem  altogether,  resulting  in  a  boundary  residual  form, 


r  B(Q)-b  \ 


(5) 


Here,  TV  is  a  (n 
that 


to)  x  n  matrix  containing  a  basis  for  the  nullspace  of  the  boundary  condition  matrix,  such 


N 


=  0. 


(6) 


In  essence,  the  Lagrange  multiplier  method  provides  a  way  to  augment  Dirichlet  boundary  conditions  with 
some  combination  of  the  equations  governing  mass,  momentum  and  energy.  The  Lagrange  multiplier  method 
is  easily  applied  to  the  node-centered  formulation  because  it  makes  use  of  the  discretized  equations  of  motion, 
R(Q ),  at  the  boundary  nodes.  As  previously  stated,  N  depends  upon  the  choice  of  Q ,  resulting  in  different 
boundary  residuals. 

A  novel  aspect  of  this  work  is  the  verification  of  boundary  conditions  through  rigorous  grid  refinement 
studies.  We  do  this  for  exact  solutions  which  are  available  for  certain  simple  configurations.  However, 
we  desire  a  more  general  method  that  is  applicable  to  a  wide  variety  of  boundary  conditions  and  flow 
equations  that  do  not  possess  exact  solutions.  To  accomplish  this,  we  extend  our  previous  work  using  the 
method  of  manufactured  solutions  (MMS)  [13],  to  include  boundary  conditions.  This  enables  the  use  of  a 
single  arbitrary  manufactured  solution  to  simultaneously  verify  the  accuracy  of  both  the  interior  scheme  and 
boundary  conditions.  This  is  accomplished  by  adding  an  MMS  source  term  to  both  the  Euler  and  boundary 
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equations: 


dQ 

dr 


+  V  •  F  =  S(x) 


B(Q)  -  b  =  Sb(x) 


For  the  Lagrange  multiplier  method  this  results  in  a  modified  boundary  residual  of  the  form 


(7) 


Rb(Q) 


B(Q)  —  b  —  Sb(x)  \ 
N(jg+R(Q)-S(xj)f- 


(8) 


The  proposed  method  of  boundary  verification  does  not  require  that  the  manufactured  solution  satisfy  the 
boundary  conditions  as  required  by  other  approaches [5].  This  greatly  simplifies  the  verification  procedure, 
and  allows  a  single  manufactured  solution  to  be  used  to  verify  any  number  of  boundary  conditions. 

We  test  three  boundary  conditions  that  make  use  of  the  Lagrange  multiplier  methodology.  These  are 
denoted  n-INVl,  nINFl,  and  n-OUTl.  Method  n-INVl  for  a  node-centered  inviscid  wall  sets  the  normal 
velocity  to  zero  and  retains  a  combination  of  the  mass,  momentum,  and  energy  equations: 


B(Q)=un,  b  =  0, 


/I  0  0  0\ 

N  =  I  0  ny  —nx  0  1 
\0  —  u  —v  1 J 


(9) 


Here,  n  =  (nx>  ny)T  is  the  outward  pointing  unit  normal  vector  at  the  surface  node,  and  un  =  nxu  +  nyv  is 
the  velocity  in  the  normal  direction. 

The  subsonic  inflow  condition  n-INFl  specifies  entropy,  s,  total  enthalpy,  h,  and  tangential  velocity, 
ut  =  —nyu  +  nxv: 


(h  \  f  hspec  \ 

S  1,  6=  Sspec  I,  N  =  (unh  utv  -  nx(h  +  \q2)  ~utu  -  ny(h  +  \q2)  un)  (10) 

J  \^t,specj 

Here,  the  subscript  spec  denotes  a  specified  value  and  q  is  the  velocity  magnitude. 

The  final  method  using  Lagrange  multipliers  is  n-OUTl  for  subsonic  outflow.  This  method  fixes  static 
pressure, 

(-u  1  0  0\ 

B{Q)=P ,  b  =  Pspea  N=  -u  0  1  0  .  (11) 

\-h  0  0  1/ 

Again,  we  emphasize  that  the  method  of  Lagrange  multipliers  provides  an  automatic  way  to  determine 
a  complete  set  of  boundary  equations.  Once  we  determine  B(Q)  and  a  set  of  variables  Q ,  the  remaining 
equations  are  fixed  through  the  nullspace,  N.  However,  the  actual  form  of  N  is  dependent  on  the  choice  of 
Q,  which  appears  to  be  arbitrary.  Further  work  is  needed  to  determine  which  set  of  Q  and  the  resulting  N 
is  best  suited  for  various  conditions. 

2.1.2  Selection  Matrix  Method 

A  more  intuitive  method  to  obtaining  a  boundary  residual  is  through  a  selection  matrix  that  picks  desired 
combinations  of  the  equations  of  motion.  Using  a  selection  matrix,  A,  we  define  the  boundary  residual  as 

Rb(Q)  =  Q(Q)  +  A  ^  +  RiQ)^.  (12) 

Here,  A  selects  certain  combinations  of  the  equations  of  motion  which  augment  other  boundary  conditions, 
H(Q).  These  conditions  may  be  Dirichlet  ( D ),  Neumann  (N),  or  extrapolated  (E): 

fl(Q)  =  Hd(Q)  +  CIn(Q)  +  flE(Q)-  (13) 
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Whereas  the  Lagrange  multiplier  approach  only  accommodates  Dirichlet  conditions,  this  framework  allows 
for  general  specification  of  virtually  any  type  of  boundary  condition  imaginable.  It  should  be  noted  that 
the  Lagrange  multiplier  method  is  a  subset  of  this  more  general  approach,  where  £Id(Q)  =  B(Q)  —  6,  and 
A  =  N  augmented  with  rows  of  zeros. 

This  more  general  form  is  also  easily  modified  for  verification  via  MMS  by  introducing  source  terms  for 
f l(Q)  and  the  equations  of  motion  themselves: 

Rb(Q)  =  fi(Q)  -  Sb(x)  +  A(^  +  R(Q)  ^  Six)')  .  (14) 


Here,  Sb(x)  acts  on  the  f2(Q)  terms,  and  S(x)  acts  on  the  governing  equations  of  motion. 

We  test  four  node-centered  boundary  conditions  that  make  use  of  the  selection  matrix  method.  These 
are  denoted  n-INV2,  n-INV3,  n-INF2  and  n-OUT2.  For  n-INV2,  a  Neumann  condition  is  set  only  for  the 
tangential  velocity  terms,  while  the  mass  and  energy  equations  are  used  directly: 


0  1 

Qd{Q)  = 

0 

Un 

,  £In(Q)  = 

dut  /  dn 

0 

'v  o  ) 

H e{Q )  =  0, 


/I  0  0  0\ 
0  0  0  0 
0  0  0  0 
\0  0  0  1/ 


(15) 


In  this  manner,  velocity  is  treated  with  a  symmetry  condition  while  mass  and  energy  are  conserved.  These 
types  of  symmetry  conditions  enforced  with  Neumann  conditions  have  been  advocated  by  many  other  re¬ 
searchers  for  slip  walls  in  inviscid  flows  [14,  15,  16]. 

Method  n-INV3  uses  Neumann  conditions  for  all  equations,  completely  omitting  any  contribution  from 
the  Euler  equations: 


(  dp/dn  \ 

(0 

0 

0 

°\ 

Vd(Q)  = 

0 

Un 

,  Hjv(Q)  = 

dut/dn 

0 

,  He(Q)  =  0,  A  = 

0 

0 

0 

0 

0 

0 

0 

0 

w 

\d(pe)/dnj 

V> 

0 

0 

0/ 

This  condition  is  strictly  a  symmetry  condition,  but  is  often  used  at  inviscid  walls  in  practice. 

For  subsonic  inflow  method  n-INF2,  the  selection  matrix  extracts  the  outgoing  characteristic  equation 
at  the  inflow  boundary,  while  three  other  conditions  (total  enthalpy,  entropy,  and  tangential  velocity)  are 
prescribed.  Here  we  wish  to  select  the  equation  associated  with  the  outgoing  un  —  c  wave.  To  accomplish 
this,  the  Euler  equations  are  first  rotated  to  a  coordinate  system  with  components  normal  (n)  and  tangential 
(t)  to  the  boundary, 


(dQ  dFx  dFy\  _  dQ'  8Fn  dFt 
\dr  +  8x  +  dy  )  dr  +  dn  +  dt 


(17) 


Here,  the  definitions  of  the  rotation  matrix,  rotated  solution,  and  flux  vectors  are 


(l  0  0  0\ 

(  P) 

(  PUn  \ 

f  Put  \ 

R  = 

0  nx  ny  0 

0  —  ny  nx  0 

^0  0  0  l) 

,  Q'  =  RQ  = 

PUn 

put 

\Pe  / 

F  — 

7  x  n  — 

PU2n+P 

pUnUt 

v  punh  y 

,  Ft  = 

putun 
pu2  +P 
\  puth  J 

Equation  17  is  then  multiplied  by  the  modal  matrix,  Mn  1 ,  containing  the  right  Eigenvectors  of  the  normal 
flux  Jacobian,  An  =  dFn/dQ' .  Using  the  fact  that  M~1AnMn  =  An  =  di&g(un,un,un  +  c,un  —  c),  leads  to 


dQ'  .  d& 

<-»  i  j- i-n  o 

or  on 


dOr 

-M~1AtMn-^  =  0, 


(19) 


where  Q'  =  Mn  1RQ.  It  is  clear  from  this  form  that  the  fourth  equation  represents  the  characteristic  wave 
equation  leaving  the  boundary.  We  wish  to  select  this  equation  as  part  of  the  boundary  specification.  This 
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can  be  accomplished  by  defining  the  selection  matrix  as 


where 


/ 


M~lR  = 


1  _  (7~1)<?2 
1  2c2 


A  = 


/0  0  0  0\ 
0  0  0  0 
0  0  0  0 
\0  0  0  1/ 

(7-1)2 


M~XR, 


(7-1)2 


Tly  'W'x 

2^2  (^g2  -  cu„)  ^s{nxc-  (7-  l)u)  ^(n^c- (7  -  l)u) 
+  CMn)  -23-(na;C+  (7-  l)u)  -  5^2  {nyC  +  (7  -  l)t>) 


7-1  \ 

c2 

0 

7-1 
2c2 
7-1  / 
2c2  ' 


(20) 


(21) 


The  matrix  M~lR  is  equivalent  to  nxM~l  +  nyM~l,  where  M'1  and  M~ 1  are  the  modal  matrices  of  the 
Cartesian  flux  Jacobians,  Ax  =  dFx/dQ  and  Ay  =  dFy/dQ,  respectively.  The  n-INF2  condition  is  completed 
by  specifying  total  enthalpy,  entropy,  and  tangential  velocity: 


Qd(Q) 


^  h  hspec  \ 

S  SSpec 
^ t,spec 

V  0  / 


^jv(Q)  =  0,  Qe{Q)  =  0. 


(22) 


The  outflow  case  n-OUT2  is  based  on  a  similar  approach  as  n-INF2,  but  selects  the  remaining  three 
characteristics  at  the  outflow,  while  fixing  the  static  pressure: 


&d(Q) 


(  °  \ 

0 

0 

\P  Pspec ) 


^n(Q)  =  0,  Qe(Q)  =  o, 


/!  0 
0  1 
0  0 
\0  0 


0 

0 

1 

0 


°\ 

0 

0 

0/ 


(23) 


2.1.3  Weak  Form  Flux  Implementation 


A  final  method  of  boundary  implementation  suitable  for  node-centered  schemes  is  through  a  modification  of 
the  flux.  Traditionally,  this  has  been  used  at  an  inviscid  wall  by  setting  the  convective  portion  to  zero  and 
retaining  only  the  pressure  terms  [2]: 


f°\ 


pnx 

pny 


(24) 


V  o 


Here,  Fn  is  the  directional  flux  at  a  boundary  face  with  normal  n  =  (nx,ny)T.  We  denote  this  method 
n-INV4  in  Table  1. 

It  is  interesting  to  note  that  this  method  does  not  explicitly  involve  a  boundary  residual,  and  thus  cannot 
be  verified  using  MMS  and  source  terms  as  the  other  methods  allow.  Nonetheless,  our  results  show  that 
this  method  performs  very  well  as  shown  through  grid  refinement  studies  using  exact  solutions.  It  is  also 
interesting  to  note  that  this  method  does  provide  a  water-tight  formulation  for  node-centered  schemes  using 
inviscid  walls. 


2.2  Cell-Centered  Boundaries 

Unlike  node-centered  schemes,  cell-centered  schemes  do  not  contain  unknowns  that  lie  directly  on  the  bound¬ 
ary.  Instead,  ghost  nodes  are  introduced,  as  shown  in  Figure  1(b).  These  ghost  nodes  do  not  lie  within  control 
volumes  and  do  not  contain  flux  balances  from  the  governing  equations.  For  this  reason,  it  is  difficult  to 
apply  any  discretized  equations  of  motion  (e.g.  mass,  momentum,  or  energy)  at  boundary  locations.  It 
is  at  these  ghost  node  locations  that  the  boundary  residual  must  be  satisfied  at  steady-state.  Because  of 
the  difficulty  of  obtaining  discretized  equations  of  motion  at  ghost  nodes,  we  limit  the  boundary  condition 
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choices  to  a  combination  of  Dirichlet,  Neumann,  or  extrapolation  conditions: 

Rb(Q)  —  H(Q)  =  ^d(Q)  +  ^n{Q)  +  £Ie(Q)- 


(25) 


These  conditions  may  also  be  verified  using  MMS  by  adding  a  source  term  to  the  boundary  residual: 


Rb{Q)  =  ^(Q)  ^  Sb(x). 


(26) 


An  important  point  is  that  for  extrapolation  conditions,  the  source  term  should  be  set  to  zero.  This  is 
because  in  the  limit  of  grid  refinement,  the  extrapolated  value  will  equal  the  manufactured  solution  exactly. 

We  test  six  cell-centered  boundary  conditions  that  make  use  of  the  ghost  node  method.  These  are  denoted 
c-INVl,  C-INV3,  c-INFl,  C-INF2,  c-OUTl,  and  C-OUT2  in  Table  1.  Method  c-INVl  for  an  inviscid  wall  sets 
the  normal  velocity  to  zero,  while  extrapolating  density,  tangential  velocity,  and  total  energy: 


nD(Q)  = 

(°) 

0 

Un 

*  Hjv(Q)  =  0,  S1e(Q)  = 

(  P  Pe  \ 
Ut  -  UtE 

0 

\pe  -  ( pe)E ) 

(27) 


Here,  the  subscript  E  refers  to  the  state  (linearly)  extrapolated  from  the  interior  of  the  domain. 

Method  C-INV3  for  an  inviscid  wall  sets  the  normal  velocity  to  zero,  but  uses  Neumann  conditions  for 
density,  tangential  velocity,  and  total  energy: 


(°\ 

/  dp/dn  \ 

^d(Q)  = 

0 

Un 

,  f In(Q)  = 

dut  /  dn 

0 

w 

yd(pe)  /  dn  J 

Me(Q)  =  0. 


(28) 


The  Neumann  conditions  may  be  more  appropriate  for  a  symmetry  condition,  but  this  is  often  used  for  an 
inviscid  wall  in  practice. 

Method  c-INFl  for  subsonic  inflow  fixes  the  thermodynamic  stagnation  state  and  incoming  tangential 
velocity,  while  extrapolating  the  velocity  normal  to  the  boundary: 


flo(Q) 


^  h  hspec  \ 


ut  -  Ut  spec  I 

V  0  / 


H]v(Q)  =  0,  &e{Q) 


(  °  \ 

0 

0 

UnE  J 


(29) 


Method  C-INF2  for  subsonic  inflow  uses  Riemann  invariants  normal  to  the  boundary,  entropy,  and  tan¬ 
gential  velocity.  The  Riemann  invariants  have  the  usual  definition  for  an  ideal  gas, 


R+  —  un  + 


2c 

7-1’ 


R  —  Un 


2  c 

7-1’ 


where  c  is  the  speed  of  sound. 


The  boundary  condition  then  becomes, 


H_d(Q) 


(  °  ^ 

So 

+ 

1 

D —  _  D  — 

J i  lxspec 

'U't  'U,t,spec 

\  &  &spec  j 

,  Hat(Q)  =  0,  £Ie(Q)  = 

o  o  o 

(30) 


(31) 


The  subsonic  outflow  condition  c-OUTl  specifies  the  static  pressure,  and  extrapolates  density  and  veloc¬ 


ity: 


Hd(Q) 


(  ° 

q  )  fijv(Q)  =  0,  £Ie(Q) 

\p  Pspec  J 

8 


(  P  —  Pe\ 
u  —  Ue 
V  -  Ve 

V  0  ) 


(32) 


The  subsonic  outflow  condition  C-OUT2  is  similar  to  C-INF2,  but  extrapolates  tangential  velocity  and 
entropy  to  be  consistent  with  the  number  of  characteristics  leaving  the  domain: 


^d(Q) 


(  o  \ 

+ 

tq  + 

JD—  _  TD  — 

''  -**•' spec 

0 

,  nN(Q)  =  o,  nE(Q)  = 

0 

ik  ~  UtE 

V  o  ) 

\  s-  SE  ) 

(33) 


The  final  condition  listed  in  Table  1  is  C-INV2.  This  condition  is  similar  to  n-INV4,  which  weakly  modifies 
the  flux  at  the  boundary  to  only  contain  the  pressure  terms,  as  in  Equation  24.  Thus  only  pressure  needs 
to  be  extrapolated.  As  in  n-INV4,  it  is  not  clear  what  the  boundary  residual  is  for  this  method.  Thus,  it  is 
difficult  to  verify  the  method  using  MMS  and  source  terms  as  the  other  methods  allow. 


3  Results 

In  this  section  we  present  results  of  the  various  boundary  condition  implementations  in  one  and  two  di¬ 
mensions  for  node-  and  cell-centered  schemes.  Four  cases  are  considered:  (1)  quasi- ID  nozzle  flow,  (2)  2D 
manufactured  Euler  solution  in  a  square  domain,  (3)  Ringleb  Flow,  and  (4)  subsonic  flow  over  a  NACA  0012 
airfoil.  For  all  cases,  we  perform  grid  refinement  studies  with  the  various  boundary  conditions  listed  in  Table 
1.  At  times,  certain  schemes  fail  to  produce  second-order  results,  in  which  cases  explanations  are  given  to 
account  for  these  observations. 


3.1  Quasi-ID  Euler  Equations 

The  quasi-ID  Euler  equations  are  a  ID  formulation  with  added  cross  sectional  area,  a( x),  useful  for  nozzle 
flows.  The  quasi-ID  Euler  equations  are  defined  as: 


dQ  dF 
dt  +  dx  ~  p 


Q  =  a 


pu 

,  F  =  a  |  pu2  +  p 
puh 


In  addition  we  assume  ideal  gases  and  constant  specific  heats. 


For  the  area,  we  use 


a(x)  =  ao  (  1  +  -cos(27ra;) 


(34) 


(35) 


(36) 


where  ao  =  2/3  to  make  the  inlet  and  outlet  area  unity.  In  addition,  the  domain  extents  are  x  €  [0, 1].  A 
2D  view  of  the  nozzle  can  be  seen  in  Figure  2. 

We  solve  the  governing  equations  using  a  node-centered  approach  to  test  inflow  and  outflow  conditions 
n-INFl,  n-INF2,  n-OUTl,  and  n-OUT2.  Reduction  of  these  methods  from  two  dimensions  to  one  dimension 
is  straightforward.  In  order  to  verify  the  accuracy  of  the  boundary  methods,  we  consider  the  exact  solution, 
which  for  fully  subsonic  flow,  may  be  found  simply  from  preserving  constant  entropy,  enthalpy,  and  mass 
flux.  For  all  tests,  M  =  0.15  was  selected  at  the  inflow  to  ensure  subsonic  flow  at  the  throat. 

While  the  quasi-ID  Euler  equations  possess  exact  solutions,  we  also  formulate  a  manufactured  solution 
to  test  our  MMS  verification  procedure  of  boundary  conditions.  The  manufactured  solution  is  chosen  as 

pMMS(x)  =  ci  +  cx\sin{ax\x) 


UMMS(x)  =  C2  +  Cx2COs(ax  2X) 
pMMS(x)  =  c3  +  cx3sin(ax3x) 


(37) 
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Figure  2:  Converging-diverging  nozzle  used  for  the  Quasi-ID  Euler  equations. 


These  constants  are  selected  such  that  the  manufactured  solution  remains  physically  meaningful  (e.g.  positive 
density).  The  period  of  the  sinusoidal  oscillation  is  set  to  approximately  20  times  the  length  scale  of  the 
problem  to  ensure  a  smooth  solution  in  the  asymptotic  range  of  convergence.  Values  for  the  constants  are 
shown  in  table  2. 


Constant 

Value 

Cl 

1.0 

Cxi 

0.15 

&X1 

0.075tt 

C2 

70.0 

Cx2 

7.0 

<1x2 

0.157T 

C3 

1.0  x  10“5 

Cx3 

2.0  x  10"4 

^x3 

0.1t r 

Table  2:  A  list  of  the  MMS  constants  used  in  Equation  37  for  the  quasi-ID  Euler  solution  verification. 


With  these  manufactured  solutions  defined,  it  is  necessary  for  us  to  determine  the  accompanying  source 
terms,  S(x)  and  Sb(x),  in  Equation  14.  The  interior  source  terms  are  computed  in  the  usual  manner  by 
setting  the  source  terms  to  the  residual  quantity  remaining  after  the  manufactured  solution  is  substituted  in 
the  governing  equations.  The  boundary  source  term  is  computed  in  a  similar  fashion,  but  considering  only 
the  boundary  equations,  such  that 

Sb  =  n(QMMS),  (38) 

where  QMMS  is  the  chosen  manufactured  solution.  For  example,  the  n-OUT2  condition,  which  specifies  the 
exit  pressure  requires  a  boundary  source  term  computed  with 


Sb  =  n(QMMS) 


(39) 


With  these  definitions  of  S  and  Sb  is  then  possible  to  solve  Equation  12  to  verify  the  interior  scheme  and 
boundary  conditions  simultaneously  using  MMS. 

The  results  of  the  grid  refinement  study  using  the  exact  and  manufactured  solutions  for  the  quasi-ID 
nozzle  are  shown  in  Figure  3.  All  methods  produce  second  order  accurate  results.  These  results  highlight 
two  important  points.  First,  the  choice  of  boundary  conditions  that  lead  to  consistent  formulations  does  not 
appear  to  be  unique.  In  this  case  the  Lagrange  multiplier  method  and  the  characteristic  matrix  selection 
method  produce  different  equations.  However,  both  equation  sets  lead  to  second-order  accuracy.  While  the 
choice  of  equations  is  not  unique,  later  we  show  examples  of  boundary  conditions  that  are  definitely  wrong 
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Figure  3:  Grid  refinement  study  for  quasi-ID  Euler  equations  using  methods  n-INFl,  n-INF2,  n-OUTl,  and 
n-OUT2  with  exact  and  manufactured  solutions. 


and  degrade  the  accuracy. 

A  second  important  point  is  that  the  MMS  procedure  provides  reliable  information  regarding  the  accuracy 
of  the  boundary  formulation.  In  this  case,  an  exact  solution  is  available.  However,  in  the  vast  majority  of 
cases  that  use  more  complex  equations  and  boundary  conditions,  exact  solutions  are  not  available.  In  these 
cases,  the  MMS  procedure  provides  a  straightforward  way  to  measure  the  accuracy  of  boundary  conditions. 
By  determining  the  proper  forms  of  S(x)  and  Sb{x),  it  is  possible  to  use  an  arbitrary  manufactured  solution 
to  verify  boundary  conditions  along  with  the  interior  scheme.  The  next  case  highlights  some  important 
considerations  when  choosing  forms  for  the  manufactured  solution. 

3.2  Manufactured  Solution  in  Square  Domain 

Here,  we  extend  the  ideas  above  to  two-dimensional  flows  and  further  explore  verification  via  manufactured 
solutions  on  a  simple  square  domain.  We  consider  the  unit  square  domain  shown  in  Figure  4.  The  domain 
is  discretized  with  perturbed  triangular  cells,  shown  in  Figure  4(a),  to  avoid  any  fortunate  cancellation  of 
solution  error  due  to  grid  regularity.  Velocity  vectors  for  the  chosen  manufactured  solution  is  shown  in  Figure 
4(b).  The  manufactured  solution  follows  a  sinusoidal  form  similar  to  the  solution  used  for  the  quasi-ID  nozzle 
flow.  An  example  for  density  is: 

p(x,  y)  =  ci  +  cxisin(axlx )  +  cylsin(ayly)  +  cxyiCOs{axyixy) .  (40) 

Here  the  constants  are  chosen  to  be  physically  meaningful  in  some  sense.  The  other  flow  variables  (u,  v,  p ) 
use  similar  forms.  For  the  computation  of  source  terms  in  multiple  dimensions  it  is  essentially  necessary  to 
use  a  symbolic  math  tool,  such  as  the  one  contained  in  Matlab. 

We  systematically  test  the  accuracy  of  various  boundary  conditions  from  Table  1  in  isolation.  To  isolate 
a  given  boundary  type,  pure  Dirichlet  conditions  were  enforced  on  all  the  other  boundaries  of  the  domain 
except  the  one  in  question.  The  inviscid  wall  cases  (INV)  were  applied  on  the  top  of  the  domain,  inflow 
(INF)  on  the  left  side,  and  outflow  (OUT)  on  the  right  side.  With  some  of  the  methods  (n-INV4  and  C-INV2) 
it  is  not  possible  to  apply  the  MMS  procedure  since  it  is  unclear  what  the  exact  governing  equations  are  for 
these  weak  forms. 

Figure  5  shows  the  results  of  the  tests  performed  for  all  the  available  MMS  boundary  types  for  node- 
and  cell-centered  methods.  Figures  5(a)  and  5(c)  show  the  results  for  the  isolated  inviscid  wall  for  both 
discretization  schemes.  The  inflow  and  outflow  results  are  combined  into  Figures  5(b)  and  5(d)  for  node-  and 
cell-centered  methods  respectively.  The  grid  refinement  study  shows  that  Dirichlet,  Neumann,  extrapolation, 
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(a)  32  X  32  perturbed  triangular  mesh 


(b)  Manufactured  solution  velocity  vectors. 


Figure  4:  Mesh  and  manufactured  solution  used  for  boundary  condition  verification. 


and  the  equations  of  motion  themselves  are  all  reliable  second-order  accurate  boundary  conditions.  The  MMS 
procedures  proves  successful  at  verifying  the  proper  implementation  of  the  interior  and  boundary  governing 
equations.  A  subtle  but  important  point  is  that  this  procedure  assumes  the  boundary  normal  vectors  are 
exact  at  the  nodes  for  the  node-centered  scheme.  This  assumption  will  be  investigated  further. 

One  important  consideration  when  using  MMS  verification  for  boundary  conditions  is  the  characteristic 
directions  of  chosen  manufactured  solution.  This  is  illustrated  by  considering  inviscid  wall  conditions.  Fig¬ 
ure  4(b)  illustrates  that  we  have  chosen  the  flow  constants  such  that  we  obtain  flow  diagonally  towards  the 
upper-right  of  the  domain.  Because  an  inviscid  wall  specifies  one  piece  of  information  (un  =  0),  this  limits 
the  location  at  which  we  can  implement  MMS  verification  for  an  inviscid  wall  to  the  top  of  the  domain.  This 
was  our  procedure  for  all  the  results  obtained  above.  Placing  the  inviscid  wall  condition  at  the  bottom  of  the 
square  domain  violates  the  wave  propagation  characteristics  of  the  manufactured  solution  because  it  would 
require  three  specified  pieces  of  information,  not  one.  The  effect  of  this  phenomenon  is  shown  in  Figure  6, 
which  shows  the  comparison  of  density  error  for  an  inviscid  wall  (c-INV3)  at  the  top  and  bottom  of  the 
domain.  It  is  noteworthy  to  mention  that  other  inviscid  wall  boundary  formulations  (n-INVl-3)  would  not 
converge  when  enforced  at  the  bottom  of  the  domain  due  to  stability  issues.  It  is  therefore  critical  that  the 
manufactured  solution  and  boundaries  are  chosen  such  that  the  characteristic  direction  is  not  violated  for 
verification  to  be  successful. 

3.3  Ringleb  Flow 

In  the  previous  section  it  was  shown  that  all  boundary  condition  methods  tested  are  second  order  accurate 
using  the  method  of  manufactured  solutions.  We  could  not  test  methods  n-INV4  and  C-INV2  because  we  do 
not  know  the  form  of  the  underlying  boundary  residual,  R},.  Here,  we  examine  the  same  methods,  including 
n-INV4  and  C-INV2,  using  Ringleb  flow  [17],  which  is  an  exact  solution  of  the  Euler  equations.  The  results 
of  this  test  case  highlight  the  importance  of  selecting  physically  correct  conditions  for  a  given  problem. 
Just  because  a  particular  boundary  condition  is  implemented  in  a  mathematically  consistent  way  (passes  an 
MMS  verification  test  for  example),  does  not  mean  that  it  is  the  correct  boundary  condition  for  a  given  set 
of  physics.  Additionally,  we  highlight  the  difficulties  of  using  derived  quantities,  such  as  entropy,  as  measures 
of  solution  error.  Finally,  this  case  also  reveals  the  importance  of  computing  boundary  normals  correctly  for 
certain  node-centered  methods. 

Ringleb  flow  describes  inviscid  compressible  flow  turning  180°  and  involves  subsonic,  transonic,  and 
supersonic  regimes.  Here,  we  focus  on  a  subsonic  portion  of  the  flow.  The  exact  solution  is  obtained  via  a 
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(c)  inviscid  wall  -  cell-centered 


(d)  inflow/outflow  -  cell-centered 


Figure  5:  Grid  refinement  study  for  a  variety  of  boundary  formulations  using  MMS. 
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Figure  6:  Example  of  the  violation  of  characteristic  directions  for  MMS  solutions 


hodograph  method  in  the  form  of  (x,y)  coordinates  parameterized  by  velocity  magnitude,  q ,  and  streamline 
constant,  k,  which  is  the  inverse  of  the  stream  function,  ^  sin#.  Here,  6  is  the  flow  angle.  The  solution 
may  be  expressed  as 


where 
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c  +  3^  +  ^^2l°g(~c 


1  7-  1  2 

1 - qz 
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(41) 

(42) 


Given  a  location  (a;,  y),  it  is  possible  to  solve  for  the  corresponding  ( q ,  fc)  pair  with  a  few  Newton  iterations. 
We  dimensionalize  the  flow  solution  with  desired  stagnation  quantities. 

A  series  of  increasingly  refined  meshes  is  used  to  verify  the  accuracy  of  boundary  conditions.  The  mesh 
and  exact  solution  for  Ringleb  flow  are  shown  in  Figure  7.  The  domain  consists  of  two  streamlines,  which 
are  treated  as  inviscid  walls,  as  well  as  an  inflow  and  outflow  at  the  bottom  and  top  portions  of  the  domain, 
as  shown  in  Figure  7(a).  We  test  each  boundary  type  in  isolation  by  enforcing  Dirichlet  conditions  on  all 
boundaries  but  the  one  in  question.  The  symmetric  density  contours  of  the  exact  solution  are  shown  in 
Figure  7(b). 

Figures  8-9  show  the  results  for  all  boundary  conditions  listed  in  Table  1.  Four  important  observations 
are  made  in  reference  to  these  figures.  First,  we  observe  in  Figures  8(a)  and  8(c)  a  substantial  decrease 
in  accuracy  for  the  n-INV2,  n-INV3  and  C-INV3  cases.  These  three  conditions,  which  were  all  verified 
as  second-order  accurate  using  MMS,  now  exhibit  first-order  accuracy  when  applied  to  Ringleb  flow.  At 
first  this  seems  surprising,  but  may  be  explained  by  the  fact  that  the  MMS  verification  only  addresses 
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(a)  16  x  32  triangular  mesh. 


(b)  Field  plot  of  Ringleb  flow  density  field. 


Figure  7:  Configuration  for  Ringleb  flow  test  case. 


the  mathematical  correctness  of  the  discretization  of  the  chosen  conditions,  not  the  appropriateness  of  the 
conditions  themselves.  In  all  three  failed  methods,  Neumann  terms  are  included  in  the  inviscid  boundary 
conditions.  The  Neumann  conditions  are  rather  arbitrary  and  do  not  correctly  represent  the  physics  for  this 
case.  In  fact,  they  are  in  direct  conflict  with  the  governing  Euler  equations,  which  explains  the  reduced 
accuracy.  This  case  illustrates  the  critical  importance  of  selecting  appropriate  boundary  conditions  for  a 
given  set  of  physics. 

A  second  observation  is  made  by  examining  Figure  8(a),  which  shows  the  effect  of  approximating  the 
boundary  normal  vector  at  nodal  locations  for  method  n-INVl.  For  Ringleb  flow,  the  exact  normal  vectors 
may  be  deduced  from  the  exact  solution,  which  provides  x  and  y  locations  along  the  inviscid  wall  streamline 
boundaries.  We  compute  these  exact  normals  and  observe  the  effect  on  the  accuracy  compared  to  the 
approximate  normal  vectors  computed  as  the  average  of  the  surrounding  face  normals  at  a  node.  We  observe 
that  using  exact  normals  and  approximate  normals  both  result  in  second-order  behavior.  However,  the  level 
of  error  using  exact  normals  is  nearly  half  the  error  using  approximate  normals.  For  other  cases  with  more 
complex  geometry,  the  effects  may  be  amplified.  Methods  to  compute  boundary  normals  consistently  for 
complex  geometry  is  an  area  of  future  work. 

A  third  observation  is  that  all  inflow  and  outflow  conditions  tested  show  sharp  second-order  accuracy, 
as  shown  in  Figures  8(b)  and  8(d).  All  of  these  methods  are  physically  meaningful  and  compatible  with 
the  current  problem.  Thus,  the  MMS  results  for  these  conditions  are  validated.  Again,  this  demonstrates 
that  the  MMS  verification  procedure  is  sufficient  if  physically  consistent  boundary  equations  are  used.  This 
is  important  because  it  provides  a  method  for  the  verification  of  more  complex  boundary  conditions  and 
interior  schemes  for  which  no  exact  solutions  are  known,  such  as  averaged  equations  with  turbulence  models, 
multiphase  flow  models,  or  combustion. 

Finally,  it  is  interesting  to  consider  the  impact  of  the  boundary  discretizations  on  other  quantities  besides 
the  conserved  variables,  such  as  entropy.  Entropy  is  sometimes  used  as  a  measure  of  error  for  inviscid  flow 
cases  for  uniform  inflow  [18,  19,  20].  Figure  9  shows  the  convergence  of  surface  entropy  for  all  inviscid 
wall  conditions  in  Table  1.  This  figure  shows  that  none  of  the  methods  are  truly  second  order  accurate 
considering  the  entropy.  Recall  that  many  of  the  methods  were  second  order  accurate  considering  the  error 
in  the  conserved  variables,  as  shown  in  Figures  8(a)  and  8(c).  This  may  be  explained  by  the  fact  that 
isentropic  flow  is  only  attained  if  all  the  conservation  laws  are  satisfied  at  all  degrees  of  freedom.  At  the 
boundary  locations,  we  do  not  solve  the  full  set  of  conservations  laws.  Thus,  entropy  may  not  demonstrate  the 
same  convergence  behavior  as  the  mesh  is  refined.  It  may  be  possible  to  formulate  boundary  conditions  which 
preserve  specific  flow  features,  such  as  isentropic  flow.  For  example,  entropy-based  boundary  conditions  have 
been  explored  by  Balakrishnan  and  Fernandez  [4] .  The  next  case  also  highlights  potential  difficulties  with 
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Figure  8:  Grid  refinement  study  for  various  boundary  formulations  using  Ringleb  flow. 
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1/h  1/h 

(a)  Inviscid  Wall  -  Node  Centered  (b)  Inviscid  Wall  -  cell  Centered 

Figure  9:  Grid  refinement  study  for  inviscid  wall  boundary  formulations  using  Ringleb  flow  using  entropy  as 
a  measure  of  error. 


using  entopy  as  a  measure  of  error. 

3.4  NACA  0012  Airfoil 

The  last  case  we  present  is  inviscid  flow  over  a  NACA  0012  airfoil  at  M  =  0.5  and  a  =  3.0°.  We  examine 
the  effects  of  the  various  inviscid  wall  implementations  in  Table  1  on  the  entropy  production  at  the  airfoil 
surface  as  well  as  other  qualitative  aspects.  We  employ  unstructured  triangular  grids  with  40,  80,  160,  and 
320  surface  points,  respectively. 

The  results  of  the  grid  refinement  test  are  shown  in  Figure  10  for  both  node-  and  cell-centered  schemes. 
For  the  node-centered  inviscid  wall  methods,  similar  trends  can  be  seen  in  this  case  as  are  seen  in  the  Ringleb 
flow  case.  Methods  n-INV2  and  n-INV3  clearly  give  first-order  results.  Initially,  methods  n-INVl  and  nlNV- 
4  appear  to  give  second-order  results.  However,  further  refinement  indicates  first-order  accuracy  in  the 
entropy,  which  is  consistent  with  the  Ringleb  flow  results.  It  is  also  clear  from  this  case,  that  symmetry-type 
conditions  are  not  suitable  for  inviscid  walls.  Qualitatively,  errors  in  the  Mach  contours  are  apparent  with 
the  symmetry-type  conditions  as  well.  Figures  ll(a)-ll(b)  show  Mach  contours  for  method  n-INVl  and 
n-INV2,  respectively.  At  the  surface,  method  n-INV2  exhibits  oddly-shaped  Mach  contours  as  a  result  of 
symmetry  enforcement. 

Like  the  node-centered  schemes,  the  entropy  produced  from  cell-centered  schemes  reduces  between  first- 
and  second-order,  as  shown  in  Figure  10(b).  However,  this  does  not  necessarily  mean  that  the  scheme  is 
first-order  in  the  conserved  variables,  as  demonstrated  bye  the  Ringleb  flow  case.  Similar  qualitative  errors  in 
the  Mach  contours  can  be  seen  for  the  node-centered  symmetry  conditions,  as  shown  in  Figures  ll(a)-ll(b). 
Similar  Errors  have  been  noted  by  Wang  [19]  and  Barth  [21].  The  extrapolation  condition  in  Figure  11(a) 
shows  smooth  Mach  contours,  while  the  symmetry  condition  in  Figure  11(b)  displays  contour  shape  errors. 


4  Conclusions 

In  this  work,  we  present  a  general  formulation  for  arbitrary  boundary  conditions  for  node-  and  cell-centered 
finite  volume  schemes.  These  conditions  may  include  any  combination  of  Dirichlet,  Neumann,  extrapolation, 
and,  in  the  case  of  node-centered  schemes,  the  equations  of  motion  themselves.  The  ease  with  which  we  may 
apply  the  conservation  equations  of  motion  is  one  advantage  of  node-centered  schemes  over  cell-centered 
schemes.  This  avoids  the  need  for  extrapolation  or  extra  Neumann  conditions  which  may  conflict  with  the 
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(b)  Cell-centered 


Figure  10:  Grid  refinement  study  for  various  boundary  formulations  using  flow  around  a  NACA  0012  airfoil 


(a)  method  n-INVl 


(b)  method  n-INV2 


(c)  method  c-INVl 


(d)  method  C-INV3 


Figure  11:  Mach  contours  for  various  inviscid  wall  treatments  of  a  NACA  0012  inviscid  airfoil  at  M  =  0.5 
a  =  3.0°. 
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interior  scheme  or  physical  configuration  of  the  problem. 

Importantly,  we  extend  the  traditional  use  of  the  method  of  manufactured  solutions  to  include  boundary 
conditions.  By  defining  a  boundary  residual  equation  to  accompany  the  interior  residual  equation,  MMS 
can  be  used  in  both  node-  and  cell-centered  contexts  to  verify  the  combined  interior-boundary  scheme.  In 
this  manner,  a  single  manufactured  solution  can  be  used  to  verify  any  number  of  boundary  conditions  con¬ 
veniently.  Computation  of  boundary  source  terms  is  straightforward,  generally  requiring  much  less  algebraic 
manipulation  than  interior  source  terms. 

A  number  of  test  cases,  including  quasi-ID  flow  and  fully  2D  flow  demonstrate  the  use  of  the  MMS 
procedure.  It  is  found  that  the  manufactured  solution  must  be  chosen  such  that  characteristic  directions 
are  respected  during  the  boundary  condition  verification  procedure.  While  the  manufactured  solution  is 
arbitrary,  boundary  locations  should  be  chosen  to  roughly  coincide  with  the  physics  of  the  manufactured 
solution.  If  this  practice  is  violated,  erroneous  results  will  be  obtained  or  non-convergence  will  result. 

While  the  correct  choice  of  boundary  conditions  for  a  given  problem  is  often  not  unique,  a  key  finding 
from  studies  involving  Ringleb  flow  and  inviscid  airfoils  is  that  certain  boundary  conditions  may  conflict  with 
the  interior  scheme  or  the  problem  configuration.  An  example  is  the  use  of  Neumann  symmetry  conditions 
for  inviscid  walls.  Even  though  the  symmetry  condition  is  verifiably  second-order  accurate  using  the  MMS 
procedure,  it  fails  to  produce  the  proper  results  when  applied  in  an  inappropriate  physical  context.  Thus,  the 
MMS  procedure  is  only  valid  when  used  with  physical  insight  to  apply  boundary  conditions  in  the  correct 
way.  If  used  correctly,  the  MMS  procedure  provides  a  method  for  the  verification  of  complex  boundary 
conditions  and  interior  schemes  for  which  no  exact  solutions  are  known.  In  addition,  we  find  it  is  best  to 
use  the  conserved  variables  directly  to  measure  the  solution  error,  as  opposed  to  derived  quantities,  such  as 
entropy.  Entropy  appears  to  converge  as  second-order  only  when  all  conservation  laws  are  satisfied,  which 
is  not  the  case  for  the  boundary  methods  investigated  in  this  paper.  It  may  be  possible  to  formulate  special 
entropy  preserving  boundary  conditions,  which  preserve  isentropic  flow. 

Future  work  will  focus  on  modifications  of  the  node-centered  approach  to  allow  for  a  water-tight  flux 
formulation.  One  possible  approach  is  to  introduce  ghost  nodes  into  the  node-centered  formulation  in  a 
manner  similar  to  the  cell-centered  formulation.  Another  important  area  of  investigation  includes  consistent 
methods  for  boundary  normal  computation.  A  possible  strategy  includes  locally  reconstructing  the  surface 
description  from  surrounding  geometry  information  using  a  least  squares  procedure.  Derivatives  of  the 
reconstructed  surface  may  then  be  used  to  estimate  normal  vectors. 
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